#--------------------------------------------------------------------------------------------------
#APPENDIX
#--------------------------------------------------------------------------------------------------
options(warn=-1)
lf<-log_open("04_empirical_app_appendix.log")
#Contour Plot:
xgrid <- seq(0, 1, by=0.05)
ygrid <- seq(-1, 1, by=0.1)
data.fit <-  expand.grid(xgrid, ygrid)

df_plot<-data.frame(data.fit, bias= data.fit$Var1*data.fit$Var2)
names(df_plot)[1] = "Probability"
names(df_plot)[2] = "Imbalance"

df_upper<-data.frame(Probability = 1-seq(0.05, 1, by=0.05),
    Imbalance_Bound = mean(df_all$Y[df_all$T==0])/seq(0.05, 1, by=0.05))
df_upper$bias = df_upper$Probability*df_upper$Imbalance_Bound


library(metR)
df_stpaul = df_all %>% filter(city == 'St. Paul')
df_stpaul$match_index = NA
match_index = 1:length(unique(block_set1[[6]]$matching_results$MatchLoopC[,1]))
df_stpaul$match_index[block_set1[[6]]$matching_results$MatchLoopC[,1]]<-match_index
df_stpaul$match_index[block_set1[[6]]$matching_results$MatchLoopC[,2]]<-match_index
df_match<-na.omit(df_stpaul)

#Observed compliance rate:
#"pg 8: The overall observed compliance rate in the experiment is 32\%.
#This implies that the first parameter should be upper bounded by $1-p=0.68$, denoted with the gray region:
df_match %>% filter(T==1) %>% summarize(mean(C==1))
non_compl_rate = as.numeric(1- df_match %>% filter(T==1) %>% summarize(mean(C==1)))

#pg 8: At this upper bound, in order for the bias to fall within the killer confounder region,
#the imbalance term across the control units would have to exceed 0.19.
results$bdim[6]/non_compl_rate


c_treatment = mean(df_match$Y[df_match$T==1 & df_match$C==1]) - mean(df_match$Y[df_match$T==1 & df_match$C==0])
p1 = df_plot %>%
        ggplot(aes(x = Probability, y = Imbalance, z = bias)) +
        annotate("rect", xmin = non_compl_rate, xmax = 1, ymin = -1, ymax = 1 , fill= "gray", alpha=0.5)+
        geom_contour_fill(breaks =c(results$bdim[6], 1000*results$bdim[6]), fill='blue', alpha=0.25)+
         geom_contour(binwidth = 0.25, col="slategray") + ggtitle("") +
         xlab("Probability of Non-Compliance in Complier Block")+ylab("Imbalance (Violation of Block PI)")+
         geom_contour(breaks=c(results$bdim[6]), col='blue', size=1)

#FIGURE A-2.1
p1 + geom_text_contour(aes(z = bias), binwidth=0.25, stroke=0.2)+ggtitle("Bias Contour Plot")+
         annotate('text', label='Observed Non-Compliance Rate (1-p=0.68)', x = 0.66, y=-0.45, angle=90)+
         geom_vline(xintercept=non_compl_rate)+
         geom_hline(yintercept=c_treatment, linetype='dashed')+
         annotate('text', label="Average Imbalance: 0.13", x=0.1, y=0.17)
ggsave('../Figures/figA2_1_contour.pdf', height=7.5, width=7.5)
#---------------------------------------------------------------------------------------------------------
#TABLE A-4.3
log_print(xtable::xtable(rel_reduc[7:12,-1], digits=1), include.rownames=FALSE)


#TABLE A-4.5
#PREP DATA FOR TABLE:
point_estimates = results %>% dplyr::select(blocking_set, cities, iv, iv_block, psw, psw_block, bdim)
se = results %>% dplyr::select(blocking_set, cities, iv_se, iv_block_se, psw_se, psw_block_se, bdim_se)
names(se) = names(point_estimates)
se$blocking_set = NA
se$cities = NA
df_estimates_all = lapply(1:nrow(point_estimates), function(x){rbind(point_estimates[x,], se[x,])}) %>% bind_rows()
df_estimates_all[,-(1:2)] = df_estimates_all[,-(1:2)]*100

#GENERATE TABLE:
log_print(xtable::xtable(df_estimates_all, digits=2), include.rownames=FALSE)

log_close()